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Abstract 

The population dynamics of a trapped Bose-Einstein condensate, subject to 
the action of an external field, is studied. This field produces a spatio-temporal 
modulation of the trapping potential with the frequency close to the transition 
frequency between the ground state and a higher energy level. For the evolution 
equations of fractional populations, a critical line is found. It is demonstrated that 
there exists a direct analogy between dynamical instability at this line and critical 
phenomena at a critical line of an averaged system. The related critical indices are 
calculated. The spatio-temporal evolution of atomic density is analyzed. 
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1 Introduction 



A dilute cloud of Bose-condensed atoms confined within a trapping potential at low tem- 
peratures can be described by a wave function satisfying the Gross-Pitaevskii equation 
[1,2]. This equation is nonlinear due to the atomic interactions through a local delta po- 
tential. The mathematical structure of the equation is that of the nonlinear Schrodinger 
equation. Stationary states of the latter, because of the confinement caused by a trapping 
potential, are restricted to discrete energy levels. These stationary eigenstates form a set 
of nonlinear modes, in analogy to linear modes that are solutions of a linear Schrodinger 
equation. The nonlinear modes can be called coherent modes, since the wave function of 
the Gross-Pitaevskii equation corresponds to a coherent state of Bose-condensed atoms. 
It is also possible to use the term topological modes, emphasizing that the wave functions 
related to different energy levels have different spatial topology. 

This modal structure of the confined condensate states is in close analogy with non- 
linear optical modes of nonlinear waveguide equations in optics [3,4]. And the methods of 
creating nonlinear coherent modes of trapped Bose atoms [5-7] are also similar to those 
employed in optics, where one uses specially prepared initial conditions or invokes an ac- 
tion of external fields. Dynamics of fractional populations, characterizing the occupation 
of coherent modes of Bose-condensed trapped atoms, displays as well many phenomena 
equivalent to those known in optics, for instance collapses, revivals, and Rabi-type oscil- 
lations. The nonlinear dynamics of processes coupling bosonic coherent modes have been 
studied in several papers considering the general case [5-7], antisymmetric modes [8,9], 
dipole topological modes in a two-component condensate [10], dark-soliton states [11], 
and vortex modes [12-19]. 

The aim of the present paper is threefold: First, we demonstrate the existence of a 
critical line for the population dynamics of Bose condensates, which is a kind of bifurca- 
tion line due to the nonlinearity of evolution equations. Second, we show that there is an 
intimate relation between the studied nonlinear dynamical system and the corresponding 
averaged system, so that the bifurcation line is an analog of the critical line, and instabil- 
ities happening at the former are the counterparts of critical phenomena occurring at the 
latter. Third, we study the spatio-temporal evolution of atomic density under the action 
of a resonant external field. 



2 Critical Phenomena 

The time-dependent Gross-Pitaevskii equation has the form 

y? , H(<p) = - ^- V 2 + U(t) + AN\ V \ 2 , (1) 

where U(r) is a trapping potential; A = Arrh 2 a s /m ; a s is a scattering length; m is mass; 
and iV is the number of particles; V is a potential of external fields. The wave function 
ip is normalized to unity, \\ip\\ = 1. The nonlinear topological modes are defined [5-7] 
as the eigenfunctions of the nonlinear Hamiltonian, that is, they are the solutions to the 
eigenproblem 

H(ip n )<Pn(r) = E n ip n (r) . 
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It is worth stressing here the principal difference between the topological coherent modes 
and collective excitations. The former are self-consistent atomic states defined by the 
nonlinear Gross-Pitaevskii equation. While the latter are the elementary excitations cor- 
responding to small deviations from a given atomic state and are described by the linear 
Bogolubov-De Gennes equations [1,2]. 

Assume that at the initial time the system was condensed to the ground-state level 
with an energy Eq. So that the initial condition to Eq. (1) is tp(r, 0) = y?o( r )- Suppose we 
wish to couple the ground state ipo with another state ipj having a higher energy Ej. The 
best way for doing this is, clearly, by switching on an external field, say V = V(r) cos cut, 
oscillating with a frequency u which is close to the transition frequency ujj = (Ej — E )/K, 
so that the detuning Alu = uo — uoj be small, \ Alu/lu\ <C I. Then one can look for a solution 
to Eq. (1) in the form 

<p(r,t) = c o (0<A)(r) e- iE ^ R + Cj(t) Vj (v) e~ iE ^ n . (2) 

The considered situation is very similar to the description of nonlinear resonant processes 
in optics [3,4]. The validity of the quasi-resonant two-level approximation (2) for the 
time-dependent Gross-Pitaevskii equation has been confirmed mathematically [5-7] and 
also proved by direct numerical simulations of the Gross-Pitaevskii equation [8,10,12-14], 
the agreement between the two- level picture and the simulations being excellent. 

The coefficients Ci(t) define the fractional level populations rii(t) = |q(£)| 2 . The 
equations for these coefficients can be obtained by substituting the presentation (2) in 
Eq. (1). It is again worth noticing the similarity of this presentation with the slowly- 
varying-amplitude approximation commonly used in optics, if one treats the factors Cj(t) 
as slow functions of time, as compared to the exponentials in Eq. (2), which implies that 
\dci/dt\ -C Ei. This approximation, being complimented by the averaging technique [20], 
permits one to slightly simplify the evolution equations for In this way [5-7], we 

come to the equations 

^ = -ia n jCo - 1 (3 c 3 e^ , ^ = -ia n Cj - \^ c e*** , (3) 
in which the transition amplitudes 

Mr 1 r 

a l3 = A - J |^(r)| 2 (2|^-(r)| 2 - |^(r)| 2 ) dv , (3 = - j V l(v)V(v)^(v) dv , 

caused by the nonlinearity and by the modulating field, respectively, are introduced, and 
the abbreviated notation a = a j, with setting a j = otj , is used. Note that the transition 
amplitude (3 is nonzero only if the potential V(r) depends on space variables. The initial 
conditions to Eqs. (3) are c (0) = 1 and Cj(0) = 0. 

Usually, one solves such evolution equations with fixed parameters a, (3 and Au>, 
keeping in mind a particular realization. Instead of this, we have studied the behavior 
of solutions to Eqs. (3) in a wide range of varying parameters. It turned out that this 
behaviour is surprisingly rich exhibiting new interesting effects. 

First of all, it is easy to notice that the number of free parameters in Eqs. (3) can be 
reduced to two by the appropriate scaling. For this purpose, we measure time in units of 
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a 1 and introduce the dimensionless parameters 




It is also evident that Eqs. (3) are invariant with respect to the inversion a — > —a, 
f3 — > —/3, Auj — > — Auj, and t — > — t. Therefore it is possible to fix the sign of one of 
the parameters, say a > 0, since the opposite case can be obtained by the inversion. For 
concreteness, we shall also keep in mind that (3 is positive. The dimensionless detuning 
is assumed to always be small, \5\ <C I. And the dimensionless transition amplitude b is 
varied in the region < b < 1. We have accomplished a careful analysis by numerically 
solving Eqs. (3). When parameter b is small, the fractional populations oscillate reminding 
the Rabi oscillations in optics, where \/3\ would play the role of the Rabi frequency. The 
amplitude of these oscillations increases with increasing b. It would be more correct to say 
that in our case there exist nonlinear Rabi oscillations, as far as Eqs. (3) differ from the 
corresponding equations for optical two-level systems by the presence of the nonlinearity 
due to interatomic interactions. This nonlinearity not only slightly modifies the Rabi-type 
oscillations of the fractional populations but, for a particular relation between parameters, 
can lead to dramatic effects. By accurately analyzing the behaviour of solutions to Eqs. 
(3), with gradually varying parameters, we have found out that there exists the critical 
line 

b + 5 ~ 0.5 , 

at which the system dynamics experiences sharp changes. This is illustrated in Figs. 1 to 
4, where the parameter b = 0.4999 is kept fixed and the critical line is crossed by varying 
the detuning 8. In Fig. 1, the detuning is zero, and the fractional populations display 
the Rabi-type oscillations. By slightly shifting the detuning to 5 — 0.0001 drastically 
changes the picture to that in Fig. 2, where the top of rij(t) and the bottom of no(t) 
become flat, and the oscillation period is approximately doubled. A tiny further variation 
of the detuning to 5 = 0.0001001 yields again drastic changes to Fig. 3, where there 
appear the upward cusps of rij(t) and the downward cusps of n (t). The following small 
increase of the detuning to 5 = 0.00011 squeezes the oscillation period twice, as is shown 
in Fig. 4. After this, making 5 larger does not result in essential qualitative changes of 
the population behaviour. All dramatic changes in dynamics occur in a tiny vicinity of 
the critical line. The same phenomena happen when crossing the line b + 5 ~ 0.5 at other 
values of parameters or if 5 is fixed but b is varied. This is demonstrated in Fig. 5 for 
varying 5 and another choice of b = 0.3. 

The unusual behaviour of the fractional populations is due to the nonlinearity of the 
evolution equations (3). Systems of nonlinear differential equations, as is known, can 
possess qualitatively different solutions for parameters differing by infinitesimally small 
values. The transfer from one type of solutions to another type, in the theory of dynamical 
systems, is, generally, termed bifurcation. At a bifurcation line, dynamical system is 
structurally unstable. 

The second aim of our paper is to show that the found instability in the considered 
dynamical system is analogous to a phase transition in a statistical system. To elucidate 
this analogy for the present case, we have to consider the time-averaged features of the 
dynamical system given by Eqs. (3). To this end, we need, first, to define an effective 
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Hamiltonian generating the evolution equations. This can be done by transforming these 
equations to the Hamiltonian form 

. dc dH e ff . dcj dH e ff 

dt del ' dt dc* 

with the effective Hamiltonian 

H eff = a n Q n + ^((3 e lA ^c*c 3 + (3* e' lA ^c*c ) . (4) 

An effective energy of the system can be defined as a time average of the effective Hamil- 
tonian (4). For this purpose, Eqs. (3) can be treated by means of the averaging technique 
[20], as it is described in detail in Ref. [5], which provides the guiding-center solutions. 
Substituting the latter in Eq. (4), together with the time-averaged fractional populations, 
results in the effective energy 

ab 2 ( b 2 \ 
Eeff= 2^[2T 2+6 ) ' 

where e is a dimensionless average frequency defined by the equation 

e\e 2 -b 2 ) = {e 2 -b 2 -e 2 5) 2 . 

The effective energy, being the time average of the effective Hamiltonian (4), characterizes 
the average features of the system. As an order parameter for this averaged system, one 
can take the difference of the time-averaged populations, 

b 2 

7] = n - rij = 1 . 

The capacity of the system to store the energy pumped in by the resonant field can be 
described by the pumping capacity Cp = dE e ff/d\/3\. The influence of the detuning 
on the order parameter is characterized by the detuning susceptibility xs — \dr]/d5\. 
Analyzing the behaviour of the introduced characteristics as functions of the parameters 
b and S, we found out that they exhibit critical phenomena at the critical line b + 5 = 0.5, 
which coincides with the bifurcation line for the dynamical system. Expanding these 
characteristics over the small relative deviation r = \b — b c \/b c from the critical point 
b c = 0.5 — 5, we obtain 

V-Vc^%i-25y/ 2 , Cf >~^-T-V\ Xs c-Lt-V 2 , (5) 

where r\ c = r)(b c ) and r — > 0. As is seen, the pumping capacity and detuning susceptibility 
display divergence at the critical point. The related critical indices for Cp, rj, and xs ar e 
equal to 1/2. These indices satisfy the known scaling relation: 

ind{Cp) + 2 indijf) + ind(xs) — 2 , 

where ind is the evident abbreviation for index. 
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In order to clarify what is the origin of the found critical effects for the studied dy- 
namical system, let us return back to Eqs. (3). Again we pass to dimensionless notation 
measuring time in units of a -1 . By means of the substitution 

c o = (^T^) ex P \ i (*> + ^t]\ , Cj = (^-^) exp ii ( qi - it] \ , (6) 



where p, q , and q\ are real functions of t, equations (3) can be reduced to the form 

dp , r dq bp . 

- = -b^l-fsmq, -=p+- 7 =^ cosq + 5, (7) 

in which q = qi — q . Note that this reduction to an autonomous dynamical system is 
valid for arbitrary detuning 5. Moreover, this system possesses the integral of motion 

I(p,Q) = \ v 2 -b^l-p 2 cosq + 5p, (8) 

which can be defined by using the initial conditions p(0) = — 1 and q(0) = corresponding 
to the conditions c (0) = 1 and Cj(0) = 0. This gives 

J(-l,0) = i-5. (9) 

The existence of the integral of motion means that the dynamical system is integrable in 
quadratures. This fact does not help much for studying the time evolution of the system, 
since the formal solutions p(t) and q(t) are expressed through rather complicated inte- 
grals, so that the system evolution, anyway, is to be analysed numerically. However, the 
property of integrability implies that the appearance of chaos in the system is impossible. 
Consequently, the observed critical effects in no way could be related to chaos. Then 
what is their origin? The answer to this question comes from the analysis of the phase 
portrait for Eqs. (7) in the rectangle defined by the inequalities —l<p<l,0<q< 2tt. 
This analysis shows that, if b + 5 < 0.5, then the motion starting at the initial point 
p(Q) = —1, q(0) = is oscillatory, with a trajectory lying always in the lower part of the 
phase rectangle, below the separatrix given by the equation 

^p 2 -byjl-p 2 cosq + 5p-b = . (10) 

When b + 5 = 0.5, the separatrix touches the initial point, so that the following motion 
occurs in the phase region above the separatrix. In this way, if we consider, under the 
given initial conditions, the parametric manifold formed by the parameters b G [0, 1] and 
(5 < 1, then the critical line b + 5 = 0.5 separates the parametric regions related to two 
different types of solutions to Eqs. (7). 



3 Spatio- Temporal Evolution 

The resonance formation of coherent topological modes can be noticed by observing the 
spatio-temporal behaviour of the atomic density. To illustrate this, we consider a cylin- 
drical trap, with the frequency ratio 

v = — (u r =uj x = uj y ) . (11) 
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Introduce dimensionless cylindrical variables 



(12) 



l r l r \ y ttiqUj- 

and the dimensionless wave function 

lpnmk(r, V, z) = /fVnmfe(r) , (13) 

in which n — 0, 1, 2, ... is the radial quantum number; m = 0, ±1, ±2, ... is the asimuthal 
quantum number; and k — 0, 1, 2, ... is the axial quantum number. The effective coupling 
parameter in dimensionless units is 

^ fra; r /^ l r ^ 

The density of particles 

p(r,f) = J%(r,f)| 2 , (15) 

according to the form (2), contains both slow functions of time, co(t) and Cj(t), as well as 
fastly oscillating exponentials. To exclude the fastly oscillating terms, we introduce the 

envelope density 

UJ / , 27r/w 

p(r,t) = — p(r,t)dt, (16) 

27T JO 

where the integration is over time explicitly entering the exponentials, while Co and Cj are 
kept fixed. Defining the dimensionless envelope density 

Z 3 

p(r,<p,z,t) = ^p(r,t), (17) 

we obtain 

p(r,(p,z,t) =n (t) \ijj (r,ip,z)\ 2 + nj(t) |^-(r, <p, z)\ 2 , (18) 

with j denoting the set {n, m, k} of quantum numbers. 

Numerical calculations for the envelope density (18) are illustrated in Figs. 6 to 
11, where it is clearly seen how strongly the spatio-temporal behaviour of the density 
(18) depends on the values of the parameters b and 5 in the vicinity of the critical line 
b+5 ~ 0.5. The behaviour of the density for the radial dipole mode (n — 1, m — 0, k — 0) 
is shown in Figs. 6 and 7. The corresponding transition frequency in units of u r , for large 
g > 1, is 

cjioo ^ 0.096(z/#) 2/5 . 

The density for the basic vortex mode (n — 0, m = 1, = 0) is given in Figs. 8 and 9. 
The related transition frequency, for g ^> 1, is 

3.424 

^ 01 ° " N)^ ' 
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And the density (18) for the axial dipole mode (n — 0, m — 0, k — 1) is shown in Figs. 
10 and 11; the transition frequency of the latter mode being 



c^ooi ~ 0.060(^) 2 / 5 . 

In all Figs. 6 to 11, we set the coupling parameter (14) g=100 and take the frequency 
ratio (11) v — 10, which corresponds to a disk-shape trap with the aspect ratio 



(19) 



as g >> 1. For the chosen values of g and v the transition frequencies are 

u ow ~ 0.216 (basic vortex mode) , 

u 001 ~ 0.951 (axial dipole mode) , 
cj 10 o — 1.521 (radial dipole mode) . 
This demonstrates that the transition frequencies are arranged in the order 

^oio < ^ooi < ^ioo , (20) 

the lowest energy level corresponding to the basic vortex mode. 

Let us note that the basic vortex mode, with the winding number \m\ = 1, sets off from 
other vortices with higher winding numbers \m\ > 2. For a given angular momentum, 
the creation of several basic vortices is more energetically profitable than the formation 
of one vortex with a higher winding number [2,16,18]. 

4 Density of Current 

When a coherent topological mode is excited in a trap, the density of current, as a function 
of space and time, also displays a peculiar behaviour. The density of current 

j(r, t) = - ^ b*(r, t) V^r, t) - <p(r, t)V<p*(r, t)} (21) 

can be rewritten as 

j(r,t) = — Im^(r,0V^(r,0. (22) 
m 

The wave function (2) is the sum 

<p(T,t)=<p {T,t) + <p j (T,t) (23) 

of the terms 

^(r,t) = Q (t)^(r)e^^, (24) 

where the index i — 0, j. 
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With the wave function (23), the density of current (21) becomes 

j(r,t) = jo(r,t) +3top(r,t) +j ir a(r,t) , (25) 
where the first term is the ground-state current density 

j (r, t) = — Im y9*(r, t) Vy? (r, f) , (26) 
m 

the second term is the topological current density 

j top (r,t) = — Im y9*(r,t)Vy^(r,t) , (27) 

due to the excited topological mode, and the third term is the interference current density 

j int (r, = ^-Im [<p* (r, t) V^(r, t) + ^(r, t) Vy9 (r, t)] , (28) 

caused by the interference between the ground-state mode and the excited topological 
mode. 

For the functions (24), one has 

y?*(r,t)V^M) = | Ci (0|V*(r)V^(r) . 
Using the substitution (6), one gets 

no(t) = |c (t)| 2 = , n 3 {t) = |c,(t)| 2 = ±±M . (29) 

Since the ground-state wave function v 5 o( r ) is rea l> 

jo(r,t) = 0. (30) 
And the topological current density (27) acquires the form 

jtop(r,t) = — nj(*)Im^(r)V^-(r) . (31) 

For example, in the case of a vortex mode, when ipj(r) ~ e tmip , with m — 0, ±1, ±2, . . ., 
one has 

Im^(r)V^-(r) = — \<Pj(r)\ 2 e v . 
Then the topological current density (31) is 

3top(T,t) = n,-(0|^(r)| 2 e v . (32) 

The total topological current is, of course, zero, 

J 3top(r,t) dr = , (33) 



provided there is no flux through the boundary. 

The interference current density (28) is a kind of the Josephson current oscillating 
with the transition frequency ouj ~ ou. If one averages over these fast oscillations, treating 
Ci{t) as slow functions of time, one gets 

7T Mnt(r,t) dt = . (34) 
Thus, the slow time variation of the summary density of current (25) is 

^ r2n/w 

27 Jo ^ dt = ' hop ^' ^ ' ^ 

being due to the topological current density (31). 



5 Discussion 

Here we have studied the resonant excitation of coherent topological modes in a trapped 
Bose-Einstein condensate. It is worth mentioning that there could be another way of 
creating those stationary modes that do not have zeros, except, may be, at infinity. 
Suppose, we wish to create a mode /(r) having no zeros. Then it is possible to perturb 
the system with the potential 

such that the given function /(r) be the ground-state wave function for the eigenproblem 

H(f) + V f (r)} /(r) = E f f(r) . 

However, this way does not allow us to form topologically different modes having different 
number of zeros, which is admissible by means of the resonant excitation. 

In conclusion, we have considered the population dynamics of a trapped Bose-Einstein 
condensate, subject to the action of a resonant spatio-temporal modulation of the trapping 
potential. A careful analysis of evolution equations has been made for the wide range of 
varying parameters. Such a variation of parameters can be easily realized by changing 
trap characteristics, varying the number and kind of condensed atoms, and by changing 
the scattering length using Feshbach resonances [1,2]. It turned out that on the manifold 
of possible parameters there exists a critical line where the evolution equations display 
structural instability. We have demonstrated that this instability for a dynamical system 
is analogous to a phase transition for a stationary averaged system. For the latter, one can 
define an order parameter, pumping capacity, and detuning susceptibility which exhibit 
critical phenomena at the critical line. The origin of this critical line is elucidated by 
showing that it divides the parametric manifold onto two regions corresponding to different 
solutions of the evolution equations. The spatio-temporal behaviour of the density of 
trapped atoms is studied for several first topological modes. 
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Figure captions 

Fig. 1. The time dependence of the fractional populations n (t) and rij(t) for b = 
0.4999 and 5 = 0. Here and in the Figures 1 to 5 the dashed line corresponds to the 
ground-state population n (t) and the solid line, to the excited- level population rij(t). 

Fig. 2. Flattening of the fractional populations, with their oscillation period being 
doubled, at b = 0.4999 and 5 = 0.0001. 

Fig. 3. The appearance of the upward cusps of rij(t) and of the downward cusps of 
n (t) for b = 0.4999 and 5 = 0.0001001. 

Fig. 4. Fractional populations versus time for b = 0.4999 and S = 0.00011. 

Fig. 5. The dynamics of fractional populations, for varying 5 and fixed b = 0.3, 
demonstrating the qualitative change of behaviour when crossing the critical line: (a) 
5 = 0.242681; (b) 5 = 0.242682; (c) 5 = 0.243000; (d) 5 = 0.3. 

Fig. 6. The spatio-temporal behaviour of the envelope density (18) as a function of 
the radial variable r, at fixed z = for the radial dipole mode (n — 1, m — 0, k — 0). 
Here and in all following figures the characteristic parameters g = 100, v = 10, and 
b = 0.4999 are fixed, while the detuning 5 is varied in the vicinity of the critical line 
b + 5 ~ 0.5. In this Figure, we set 5 = 0.00010. Time is measured in dimensional units 
explained in the text: t — (solid line); t = 3 (long-dashed line); t — 15 (short-dashed 
line) . 

Fig. 7. The same as in Fig. 6, but for the detuning 5 = 0.00011 and for the moments 
of time: t — (solid line); t — 10 (long-dashed line); t — 26 (short-dashed line). At the 
time t — 26 the system is in the pure radial dipole state. 

Fig. 8. Excitation of the basic vortex mode (n — 0, m — 1, k — 0). The envelope 
density is shown as a function of the radial variable r for z = 0. The detuning is S — 
0.00010. The moments of time are: t — (solid line); t — 3 (long-dashed line); t — 15 
(short-dashed line). 

Fig. 9. The same as in Fig. 8, but for the detuning 5 = 0.00011. The moments of 
time are: t — (solid line); t — 10 (long-dashed line); t — 26 (short-dashed line). At 
t — 26 the system is in the pure vortex state. 

Fig. 10. Excitation of the axial dipole mode (n = 0, m — 0, k — 1). The envelope 
density as a function of the axial variable z at the point r = 0. Here 5 = 0.00010. The 
time moments are: t — (solid line); t — 3 (long-dashed line); t = 15 (short-dashed line). 

Fig. 11. The same as in Fig. 10, but for the detuning 5 = 0.00011 and for the times: 
t — (solid line); t — 20 (long-dashed line); t = 26 (short-dashed line). At the moment 
t = 26 the system is in the pure axial dipole state. 
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